Expansion Around the Mean-Field Solution of the Bak-Sneppen Model 
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We study a recently proposed equation for the avalanche distribution in the Bak-Sneppen model. 
We demonstrate that this equation indirectly relates t, the exponent for the power law distribution 
of avalanche sizes, to D, the fractal dimension of an avalanche cluster. We compute this relation 
numerically and approximate it analytically up to the second order of expansion around the mean 
field exponents. Our results are consistent with Monte Carlo simulations of Bak-Sneppen model in 
one and two dimensions. 
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The Bak-Sneppen (BS) model jl| has become one of 
the paradigms of Self-Organized Criticality (SOC) @. 
The rules of its dynamics are very simple: the state of 
the model is completely defined by L d numbers fi ar- 
ranged on a (i-dimensional lattice of size L. At every 
time step the smallest of these numbers and its 2d near- 
est neighbors are replaced with new uncorrelated ran- 
dom numbers, drawn from some distribution V(f). This 
"minimalistic" dynamics results in a remarkably rich and 
interesting behavior. In fact, there exists a whole class of 
models called extremal models Q , which evolve according 
to similar rules, and share many similar features with the 
BS model. In all these models the update happens only 
at the site carrying the global minimum of some vari- 
able. The oldest, and perhaps the most widely known of 
these models is Invasion Percolation The BS model, 
being the simplest and the most analytically treatable 
extremal model, occupies the place of an "Ising model" 
in this class. 

The self-organized critical nature of the BS model (as 
well as of other extremal models) is revealed in its ability 
to naturally evolve towards a stationary state where al- 
most all the variables fi are above a critical threshold f c . 
The dynamics in the stationary state is characterized by 
scale-free bursts of activity or avalanches, which form a 
hierarchical structure of sub-avalanches within big- 
ger avalanches. The introduction of an auxiliary param- 
eter / [|| allows one to describe the system within the 
paradigms of standard critical phenomena. Indeed the 
distribution P(s, /) of avalanche sizes s close to f c , has 
the same qualitative behavior of the cluster distribution 
of percolation || above and below the critical thresh- 
old p c : For / < f c , P(s,f) has a finite cut-off, remi- 
niscent of an undercritical system. As f — » f c the cut- 
off diverges and a scale-free distribution P(s, f c ) ~ s~ T 
emerges. In the overcritical regime / > f c there is a 
non-zero probability to start an infinite avalanche, but 
all finite avalanches again have a finite cut-off. Scaling 
arguments || allow one to derive all critical exponents 
of a general extremal model in terms of only two inde- 
pendent ones, say r and D - the fractal dimension of the 



avalanche cluster. 

In order to compute these two remaining exponents one 
has to resort to methods which go beyond scaling argu- 
ments. Apart from the solution of the mean field case [Q, 
and a real space renormalization group approach Q for 
d = 1, a systematic theory to compute the BS exponents 
is still lacking. A promising step in this direction was re- 
cently taken by one of us |9| with the introduction of an 
exact equation for the avalanche distribution. Hereafter 
we will refer to this equation as the Avalanche Hierarchy 
Equation (AHE). It was shown that inside the AHE is 
hidden an infinite series of equations, relating different 
moments of the avalanche size distribution. 

In this letter we demonstrate that, as it was conjec- 
tured in H, the AHE indirectly relates these two expo- 
nents, thus reducing the number of independent critical 
exponents in the BS model to just one. Contrary to sim- 
ple rational relations based on scaling arguments [|| , this 
exponent relation is highly non-trivial. First we display 
the numerical solution of the AHE. Then we perform a 
perturbative "e" -expansion around the mean field solu- 
tion, up to the second order in e. The numerical solu- 
tion of AHE is in agreement with the MC simulations in 
d = 1, 2 p,^0| and is well approximated by the results of 
the e expansion up to second order. They constitute a 
significant step forward towards the full solution of the 
BS model. However, the unusual type of the expansion 
around the mean-field exponents leaves open the ques- 
tion of the upper critical dimension d c in the model. It 
also does not answer the question about the geometrical, 
fractal properties of the avalanche cluster. Instead, given 
an avalanche fractal dimension D it enables one to de- 
rive the power r of the avalanche distribution. Similarly, 
in ordinary percolation the power of cluster distribution 
r is related to the cluster's fractal dimension D via a 
hyperscaling relation r = 1 + d/D [|[. 

Following ref. ||, let us consider the exponential dis- 
tribution V(f) = e~~f , f > 0. This simplifies the expres- 
sions without loss of generality |llj . To define avalanches 
one records the signal of the model, i.e. the value of the 
global minimal number f m in{t) as a function of time t. 
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Then for every value of an auxiliary parameter /, an /- 
avalanche of size (temporal duration) s is defined as a 
sequence of s — 1 successive events, when f m i n (t) < /, 
confined between two events, when f m in(t) > /• ln 
other words, the events, when f m in(t) > /, divide the 
time axis into a series of avalanches, following one an- 
other. The AHE is an equation for the probability dis- 
tribution P(s, /) of /-avalanche sizes s, and it reads 
@: d f P( S J) = E S s ;liR d (si)P(suf)P(s - sij) - 
R d (s)P(s, /). Here R d (s) is the average number of dis- 
tinct sites updated at least once during an avalanche 
of size s. The equation describes how the sequence of 
avalanches changes when the value of / is raised by an 
infinitesimal amount df. The first term on the RHS 
describes the gain of avalanches of size s due to merg- 
ing of two consecutive avalanches of size s± and s — s±. 
Such merging occurs when the value of f m i n (t), termi- 
nating the first avalanche, happens to be in the interval 
[/, / + df]. The factor R d (si) comes from the fact that 
when avalanches merge, the active site starting the sec- 
ond avalanche must be one of the sites updated in the 
first avalanche. The bigger is the region covered by an 
avalanche, the better are its chances to merge with the 
one directly following it. The second term in the RHS 
is the "loss term" due to avalanches of size s merging to 
form a larger avalanche. 

To proceed further one needs to introduce the scaling 
ansatz R d (s) ~ for the number of updated sites. This 
power law relation is a consequence of the spatio-temporal 
fractal structure of the avalanches in the BS model |fL2| . 
The exponent fi(d), which is an independent "input" vari- 
able in the AHE, depends on the dimensionality of the 
model and the fractal structure of the avalanche. Phys- 
ically, the exponent fi relates the volume of the spatial 
projection of an avalanche cluster to its temporal dura- 
tion s. If this spatial projection is a dense object with a 
fractal dimension equal to the dimension of space d, fi is 
given by dj D. In this expression D is the fractal dimen- 
sion of the avalanche || defined through s = R D . This is 
known to be true in d = 1, where the connected nature 
of an avalanche cluster ensures the compactness (absence 
of holes) of its projection. In d — 2 the projection of 
the avalanche can have holes, but still it was numerically 
found to be dense (i.e have a fractal dimension d) ||] . It is 
clear that, as the dimensionality of space is increased, the 
exponent fi should approach 1, since multiple updates of 
the same site become less and less likely and the volume 
of the projection should be closer and closer to the total 
volume (2d+ l)s of the avalanche itself. The "hyperscal- 
ing" relation fi — d/D will be clearly violated for d > d c , 
where one has D = 4, fi = 1 |p|,|l3|. From this it follows 
that d c > 4 for the BS model. 

The introduction of the "phcnomcnological" exponent 
fi closes the AHE, which then reads 



s-l 

dfP( S , f)=Y, slP( Sl ,f)P( S - 81, f) ~ S»P(S, /). (1) 

81 = 1 

The solution of equation (Q) exhibits a power law behav- 
ior P(s, f) ~ s~ T when / is at its critical value f c . Close 
to f c it takes a scaling form 

P(s,f) = s- T F(s°Af), (2) 

where A/ = f c — f. The exponents r, a, and fi are 
related through t = 1 + [i — <x |J. Perhaps a more fa- 
miliar form of this exponent relation involves the cor- 
relation length exponent v = 1/aD. The relation then 
becomes t = l + (d—l/v)/D. It has been conjectured || 
that equation (|l]) also indirectly relates the exponents fi 
and r. In order to check this conjecture we numerically 
integrated Eq. (|l|) forward in / with the initial condi- 
tion P(s, 0) = S St i for several values of fi In order 
to locate the critical point / c (/■*)> a least square fit of 
logP(s, /) vs logs was performed runtime for each value 
of /. The value X 2 (f) °f the sum of the squared dis- 
tances from the fit drops almost to zero in a very narrow 
region (see Fig. [j]), which then allows for a very precise 
estimate of f c (fi) and r(/x). The results for the latter 
are shown in figure (□). These results are in a perfect 
agreement with Monte Carlo estimates of r in one and 
two dimensions |Q. In d = 1, fi = 1/D = 0.411(2) and 
t = 1.07(1), while the results of numerical integration of 
(0) give t(0.4) = 1.058. In d = 2 the Monte Carlo results 
fi = 2/D = 0.685(5), r = 1.245(10) are also consistent 
with our relation giving t(0.7) = 1.238. This confirms 
that Eq. (Q) indeed contains a novel non-trivial relation 
between r and fi. 

In order to address this relation analytically, let us 
take the Laplace transform of Eq.(|l|) (§. The AHE, with 
P(<*, f) = T,T=i p ( s > /)e" QS , reads 

OO 

d f ln[l - p(a, /)] = ]T P(s, f)s»e- as . (3) 

s=l 

p(a, f) has the scaling form given by: 

p(a,f) = l-a T - 1 h(Af/a°). (4) 

This scaling form follows from Eq.(||) and the scal- 
ing functions are related through h{x) = f™[F(0) - 
F{xy°)e-y]y- T dy. 

The scaling function h(x) (as well as F(x)) is analytic 
at x = 0. Its large \x\ asymptotics is determined by the 
fact that at any A/ ^ 0, p(a, f) is analytic in a, since 
P (aJ) = 1 - PooCf) + (s) f a + (s 2 ) f a 2 + . . ., and all 
the moments of P(s, f) are finite except at the critical 
point. Here Poo(f) is the probability to start an infinite 
avalanche, and Poo(/) = for / < f c . Matching the 
expected behavior of p(u, f) to its scaling form one gets: 

oo 

h(±\x\) = \x\^-°V°Y, b t\x\- k/ ' T - (5) 

fc=0 
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In other words, h{x)\x\(' y ~^l' J for \x\ 1, must be an 
analytic function of |a;| -1 /' 7 . The coefficients for k > 
are related to the amplitudes of the diverging moments 
through (s k ) f = (-l) k+1 bf\f c - i n the 

under-critical (b£ , / < f c ) and in the over-critical regime 
(b7 , f > f c )- As we will see later, it is the condition that 
h(x) has the desired asymptotic form (ph, which fixes the 
value of a for a given /x. 

With the scaling ansatz (g), the LHS of Eq.(g) be- 
comes -a- <r h'(x)/h(x), where as = A//a CT . The RHS 
needs some more work: using s ll e~ as — - 



and the identity s M 



f 



*dt/r(l — /i), one can 



express the RHS of Eq. @ in terms of an integral in- 
volving p(a, /). Matching the powers of a in the result- 
ing equation for h(x) one gets once more the well known 
exponent relation: 



r = 1 + fi — a. 
After eliminating r, equation (||) finally reads: 



(6) 



h'(x) = 



h{x) 

r(i - m) y 



1 xzh'(xz) 
dz- 



-h(xz) 



(I - z^ivy 



We were not able to solve equation (g) and find the exact 
relation <r(/i). However, we can explicitly solve AHE for 
fj, = 1. This corresponds to the mean field version of the 
BS model, which has been studied in detail 0. We will 
rederive their results using our approach. Our strategy 
will then be to perform a systematic e - expansion around 
the mean-field solution, where e = I — /i. This clearly 
differs from the standard e = d — d c expansion (note that 
the upper critical dimension d c for the BS model is still 
an open issue), since the dimensionality of the system 
does not enter directly into our discussion. 

As jj, — > 1, the integral in Eq. (0) diverges for z ~ 
1, but so does T(l — /i). This implies that the factor 
[oT(l - /i)(I - z 1 '" 7 )^]- 1 behaves as a 5{z - 1) function. 
For fi = 1 Eq. (0) reduces to 



h'(x) = axh(x)h'(x) - (I - a)h 2 (x). 



(8) 



Its solution reads - 
integration constant. Eq. 
h(x) ~ aT^C + L>:r- ^ - 
lution of Eq.® only if C 



= <Zoj with ao an 
^) implies a large x behavior 
. . .), compatible with the so- 
= 1 and a = I - a = 1/2. For 



a 



1 one recovers a mean field solution [Q: 



V4 + x 2 



(9) 



The above derivation demonstrates that the knowledge 
of the whole scaling function h{x) is necessary in order 
to find the exponent a. 

To proceed beyond the mean field case we perform a 
1— [i = e - expansion of Eq. (Q) around the /i = 1 solution 
(|). We put a = -i +ce, where c is to be determined later. 
It is convenient to change variables to z = h^(x) and 



to set h(x) = z[l + e<f>(z) + 0(e 2 )]. Keeping only terms 
linear in e one gets an equation for <fi: 

iz(l + z 2 )8 z <P = cj>- I - /2) + (1 + 2c)z 2 + 2 In 

where ip( x ) is the logarithmic derivative of T(x). This 
equation has to be solved with the boundary condition 
(f>(z = 1) = (i.e. h(0) — I). After some algebra one 
gets: 



(z) = A 



1-3? 
1 + z 2 



2 In 



1 + z 2 2 + (2 - 4c)z s 



1 



Inz (10) 



where A = 2 + ^>(l/2) = 0.03649 .... The value of c is set 
by the requirement that 4>(z) must give rise to the desired 
asymptotic behavior of h(x). The singular behavior at 
z ~ 1/x — > must be matched to the asymptotics of 
h(x) for x — > oo. To order e, Eq. (||) requires that 



[to order e°, /(y) = (yT 



5 /(a; 2 + 4ce ) ; where /(y) is analytic at y = 



3y — f)/2]. Expanding this 
relation to order e one finds that the singular part of 
(7) 4>{z) must have exactly the form 2[1 + 2cz 2 /(l + z 2 )] Inz. 
The only value of c which matches this requirement to 
the last term in the RHS of Eq. @ is c = 0. Note that 
the whole asymptotic behavior, and not just its leading 
part, is necessary to determine c. 

This concludes the first order of the expansion in e. 
We have found that in the first order in e the critical 
exponent a did not change. The exponent relation (|^) 
then gives r = 3/2 — e. Finally, the analytic form of the 
scaling function h(x), containing all information about 
the amplitudes of avalanche moments, is given by 



h(x) 



V4~- 



2i-e r 

1 + T 



f 



i Ax 



+ 0(e 2 ). 



The extension of this procedure to higher orders in e 
is straightforward, even though it involves much heavier 
algebra. Skipping the details EB| , up to second order in 
e we find 



1 4 

2-3(7 + 1*2- 



0.5 
1.5 



0.3605e 2 



0(e 3 ) 



0(e 3 



e + 0.3605e 2 + O(e 3 ). 



(11) 
(12) 



Here 7 ~ 0.5772 is the Euler's constant. The explicit 
expression for h{x) at this order is not particularly il- 
luminating, so we refrain to display it here. As seen in 
Fig. ^|, the expansion up to the first two orders is in ex- 
cellent agreement with numerical data down to p, ~ 0.6 
(c«0.4). 

On the other side, Fig. || seems to suggest a singular 
behavior of t(/z) as /1 — > 0. The specialty of fi — can 
be understood by observing that in this case P(s, /) does 
not obey scaling. Indeed, fi = corresponds to a trivial 
model with only one constantly updated site, which can 
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be considered as a O-dimensional lattice. The probabil- 
ity of /-avalanches of size s is trivially derived from the 
probability P(f m i n {t) < /) = 1 — that the signal is 
below /: P{s, /) = e~ f (1 - e^)* -1 - This is indeed the 
solution of Eq. (Q) with s M = 1 . There is no phase transi- 
tion (numerically we found / c (/i) ~ 1/fx — ► oo as /i — * 0) 
and the avalanche distribution always has an exponen- 
tial cutoff. This suggests that d — can be interpreted 
as the lower critical dimension for the BS model (note 
that P(s, /) in the d = BS model is very similar to the 
cluster size distribution in the d = 1 percolation [pi). 

In conclusion, we have shown that the avalanche hier- 
archy equation introduced by one of us in [H yields a new 
relation between exponents in the Bak-Sneppen model, 
thus reducing the number of independent exponents to 
just one. This relation expresses r, the power law expo- 
nent in the avalanche probability distribution, in terms of 
D, the mass dimension of an avalanche cluster. We were 
able to perform a systematic expansion of this relation 
around the mean field exponents, carried to the second 
order in this work. The success of this approach suggests 
that a complete e = d — d c expansion for the BS model 
could be possible. The accomplishment of this task, how- 
ever, calls for a systematic study of the BS model in high 
dimensions and for identification of the upper critical di- 
mensionality d c . 
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FIG. 1. The plot of X 2 and r vs / for fi = 0.411 and s max = 521, 1024 and 2048. 
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FIG. 2. The relation r(/i): numerical solution of Eq.(l) (□, with s max = 1024), expansion up to order 1 — /i (dashed line) 
and (1 — /i) 2 (full line). The results of the Monte Carlo numerical simulations in d = 1, 2 [3] and the mean field result [7] are 
also shown. 



5 



